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SUMMARY 

This  paper  presents  a  composite  muitigrid  method  and  its  application  to  a 
geometrically  complex  flow.  The  treatment  of  the  interior  boundary  boundary 
conditions  within  a  composite  multigrid  strategy  is  described  in  detail  for  a  1-D 
model  equation.  For  the  Navier-Stokes  equations,  a  staggered  grid  technique  is 
adopted  for  spatial  discretization  and  a  fractional  step  method  is  used  for  the 
time  advance.  Lid  driven  cavity  flows  are  used  to  demonstrate  the  effectiveness 
of  the  method. 

KEY  WORDS 

Unsteady  Navier-Stokes,  composite  multigrid  method,  fractional  step  method, 
staggered  grid. 

1.  INTRODUCTION 

Many  methods  for  numerical  simulation  of  fluid  flows  have  been  proposed 
and  CFD  (Computational  Fluid  Dynamics)  has  become  a  practical  design  tool. 
However,  problems  remain  with  regard  to  computational  efficiency,  accuracy, 
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and  turbulence  modeling.  The  most  difficult  problem  may  be  mesh  generation 
for  geometrically  complicated  domains. 

Many  schemes  have  been  devised  to  cope  with  these  issues.  The  multigrid 
method  [1]  is  one  of  the  most  efficient  schemes  for  elliptic  problems  and  has 
been  applied  to  the  Navier-Stokes  equations.  To  deal  with  geometric  complex¬ 
ity,  domain  decomposition  has  been  proposed:  its  origins  go  back  to  Schwarz 
[2],  The  idea  is  that  a  complicated  domain  be  decomposed  into  subdomains, 
whose  geometry  is  simple  enough  to  be  easily  gridded. 

In  the  present  work,  we  apply  a  combination  of  these  two  methods,  that  is,  a. 
composite  muitigrid  strategy,  to  flows  in  geometrically  complicated  domains.  In 
particular,  we  treat  2-D  unsteady  laminar  flows.  The  Navier-Stokes  equations 
are  discretized  using  second  order  central  differencing  on  a  staggered  grid  in 
space  and  a  fractional  step  time  advance  method.  The  velocity  components  are 
advanced  explicitly  and  the  pressure  is  obtained  by  solving  a  Poisson  equation 
using  a  composite  multigrid  method.  The  momentum  equations  are  integrated 
independently  on  each  subgrid.  Interpolation  on  the  composite  grid  is  accom¬ 
plished  with  a  Coons  patch  method  [3]. 

Prior  to  solving  the  Navier-Stokes  equations,  we  investigate  the  effectiveness 
of  the  composite  multigrid  technique  for  a  1-D  model  equation.  We  also  discuss 
the  interior  boundary  conditions. 

After  demonstrating  the  properties  of  the  composite  multigrid  method,  we 
appply  it  to  the  Navier-Stokes  equations.  First,  we  check  the  accuracy  of  the 
method  for  lid-driven  cavity  flow  at  Re  =  3200  and  evaluate  the  error.  Then, 
we  simulate  geometrically  complex  cavity  flows. 

2.  THE  GOVERNING  EQUATIONS 


The  2-D  unsteady  incompressible  Navier-Stokes  (N-S)  equations  can  be  writ¬ 
ten  in  the  following  form  in  a  curvilinear  coordinate  system. 


dq  d£  d£  dHj_  £//,  dP  8Q 
dt  +  +  dr\  d£  +  dq  dq 

where 

9  =  9  =  [«*  »]T>  J  -  £rb»  -  iyHt 


(1) 

(2) 


and  u  and  v  are  the  Cartesian  velocity  components  along  the  z  and  y  axes  and 
(  and  q  are  arbitrary  curvilinear  coordinates.  J  is  the  Jacobian.  Subscripts  x 
and  y  stand  for  derivatives  with  respect  to  x  and  y. 

The  second  and  third  terms  of  LHS  (1)  are  the  convection  terms  so  that  E  and 
F  are 


E  =  Uq,  F  =  Vq 


(3a) 


where 


V  = 


V  =  UTfc  +  VT)y . 


(36) 


arc  the  contravariant  velocity  components;  however,  we  shall  use  the  Cartesian 
velocity  components  as  the  primary  variables. 

The  fourth  and  fifth  terms  are  the  pressure  gradient  terms;  H(  and  Hn  can  be 
written  as: 


where  p  is  the  pressure. 

The  last  two  terms  are  the  viscous  terms;  P  and  Q  have  the  form. 


p  _  1  / (tl  +  £?)««  +  (£*»?*  +  £y»fyK  A 

J  Re  \  ((r  +  +  ((s'!*  +  ZyVy)Vn  ) 

n  _  1  ( ( ’ll  +  +  (6>Jr  +£*ifc)u«  \  /  =  > 

V  J  Re  V  (*?*  +  VI'i  +  1 

where  lie  is  the  Reynolds  number. 

The  above  equations  have  to  be  solved  simultanously  with  the  continuity  equa¬ 
tion: 


3.  NUMERICAL  SCHEME 


3.1  Differencing  in  space  and  time 

The  momentum  equations  for  the  velocity  components  u  and  v  can  be  written 
j  +Lu(u,v,Uf,  v„,p)  =  0  (7) 

)  +  Lv(u,v,u(,un,v(,vt„p)  =  0  (8) 

where  Lu  and  V  are  the  sum  of  the  convection,  pressure  gradient,  and  viscous 
diffusion  operators. 

Multiplying  Eq  (7)  by  £x  or  rjx  and  multiplying  Eq  (8)  by  £y  or  rjy  and  summing, 
we  get  the  equations  for  U  and  V . 

Jh?  00  +  (£)  (£)  = 0  <9) 

To  discretize  this  equation  in  space,  we  use  second  order  central  differencing 
on  a  staggered  grid  in  which  the  pressure  node  is  located  at  the  cell  center 
and  the  contravariant  velocity  components  U  and  V  are  located  on  the  cell 


d_/u 

dt\J 

l(l 

dt\J 
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boundaries  as  shown  in  Fig.  1.  Control  volume  I  is  used  for  the  x-momentum 
equation  and  control  volume  II  for  the  y-momentum  equation.  In  order  to 
evaluate  ti  and  v  at  the  velocity  nodes,  we  must  evaluate  U  and  V  at  those 
nodes,  but,  on  a  staggered  grid,  U  and  V  are  known  at  different  points.  Maliska 
and  Raithby  [4]  used  the  average  of  the  values  at  four  surrounding  points  to 
obtain  the  velocity  component  on  a  node  of  the  other  component.  For  example, 
to  evaluate  V  at  the  point  A  of  Fig.  1,  they  took  an  average  of  V  at  B,C,D 
and  E.  We  adopt  this  method. 

After  evaluating  U  and  V,  we  get  u  and  v  from  the  following  relation. 

(»Mt  sroo 

For  discretization  of  the  convection  terms,  we  use  the  QUICK  scheme  [5]. 

Time  advance  is  carried  out  using  the  fractional  step  method.  Let  un,  vn 
and  pn  be  the  velocity  and  pressure  at  time  step  n  and  assume  that  u"  and  r" 
satisfy  the  continuity  equation. 

First,  we  evaluate  Un+1  and  Vn+l  by  explicit  integration  of  Eq  (9). 


J 


+Lv(un,vn . Pn) 


This  is  the  predictor  of  the  MacCormack  scheme  [6]  and  many  others. 
The  corrector  for  the  MacCormack  scheme  is 


(ID 


(12) 


. ^(i) 

. pn)(Jy))]  0 

where  un+1,»n+1  are  computed  from  Un+1,Vn+l  using  Eq  (10). 

Using  Eq  (11),  Eq  (12)  can  be  rewritten  as 

. ■<>(£) 

+iV,r,‘ . p")(^)+£>::Tr,t>::?7,..,pB)(j') 

+r(u^,^,...,p")^)]  (13) 

The  new  velocity  field  u,  v  does  not  satisfy  the  continuity  equation  ( 10)  so  we 
introduce  a  pressure  correction  and  compute  the  new  velocity  field  £/n+1,  Vn+1 
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from 


fun+l\ 

1  ( 

l  vn+l J 

1  tUJ-  l 

^(aff"))  (%) 


(14) 


where  AH(U  =  ■ff”,]*'1  -  #£,,  etc.  are  determined  by  forcing  {/"+1  and  Vn+1 
to  satisfy  the  continuity  equation.  In  this  way,  we  obtain  a  pressure  equation 
(strictly  speaking,  a  pressure  increment  equation). 


+!  [”'  + 1 ;<*'«•>) + ”•  + 

-*($)+*($) 

This  is  a  Poisson  equation  for  A p  =  pn+1  —  pn  and  can  be  solved  by  the  multi¬ 
grid  method.  Because  we  use  a  staggered  grid,  no  boundary  condition  is  needed 
for  the  pressure  equation  provided  the  coordinate  system  is  orthogonal  at  the 
boundary;  this  is  the  case  in  the  present  work.  However,  the  discrete  system  of 
equations  is  singular.  For  a  simplicity  of  coding,  we  use  a  fictitious  pressure  node 
outside  the  computational  domain  to  eliminate  the  singularity.  Furthermore,  to 
ensure  uniqueness,  we  set  the  solution  at  an  arbitrary  point  (for  instance,  at 
(4,  J])  =  (1/2, 1/2))  to  zero.  Other  methods  of  de-si ngularizing  the  system  are 
available. 


3.2  Multigrid  method  for  the  pressure 

Multigrid  [1]  is  a  well  known  method  for  solving  elliptic  equations.  We  apply 
it  only  to  the  pressure  equation.  The  essential  idea  of  the  multigrid  method  is 
reduction  of  the  error  using  several  grids  of  different  sizes;  high  frequency  com¬ 
ponents  of  the  error  are  removed  by  smoothing  on  a  fine  grid  and  low  frequency 
components  are  damped  on  coarser  grids. 

The  system  of  linear  equations  (15)  can  be  written 

L4>  =  f  (16) 

where  is  the  solution  and  /  is  the  source  (forcing)  term. 

After  we  relax  Eq  (16)  using  suitable  smoother  on  the  finest  grid,  we  get  an 
approximate  solution  4>k  and 


Lkjk  =  f-Rk 


(17) 
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where  R  is  the  residual  and  superscript  k  stands  for  the  grid  level;  k  =  kmax 
corresponds  to  the  finest  grid  and  k  —  1,  a  grid  twice  as  coarse.  The  smoother 
should  reduce  the  high  frequency  components  of  the  error  rapidly.  Subtracting 
Eq  (17)  from  Eq  (16),  we  get  an  equation  for  the  error  ek : 

Lkik  =Rk  ,  ek  =  4>  -  4>k  (18) 

Next,  we  remove  lower  frequency  components  of  e.  To  do  that,  Eq  (18)  is 
solved  on  grid  k  —  l. 

Lk-1ek~l  =  Rk~'  ,  Rk -1  -  I^~1Rk  (19) 

where  Ik~l  is  a  restriction  operator. 

After  smoothing  e  on  the  coarse  grid,  e  is  interpolated  onto  the  fine  grid  by 
interpolation  or  prolongation: 


e*  =  iLi**-1  (20) 

where  Jk_ ,  is  an  interpolation  operator. 

The  solution  is  then  updated: 

4>k  =  <t>ou  +  tk  (21). 

In  the  present  work,  since  we  use  a  staggered  grid  system,  /£-1  and  1  are 
defined  as  follows. 

=  (tifjj  +  tij+ljj  +  +  ^*/+l.j/+l)/4  (22) 

where  if  =  2 ic  —  1  and  jf  —  2 jc  —  1 

Ik-itiut  =  +  3  ^e+1Je  +  3<ie+l  +  ^c+Uc+1)/16  (23) 

where  ic  =  1  +  if/2  and  jc  =  1  +  jf /2  (see  Fig.  2). 

The  method  is  readily  extended  to  larger  numbers  of  grids. 

4.  DOMAIN  DECOMPOSITION  TECHNIQUE 

In  order  to  solve  the  Navier-Stokes  equations,  we  have  to  divide  the  domain 
using  a  suitable  grid  system.  For  flow  fields  that  do  not  have  simple  geometry, 
covering  the  entire  domain  with  a  single  grid  is  difficult. 

Domain  decomposition  is  a  method  of  coping  with  this  problem.  In  this 
technique,  a  geometrically  complicated  domain  is  divided  into  several  simpler 
ones.  Schwarz  [2]  proposed  an  alternating  solution  procedure  for  elliptic  prob¬ 
lems.  Since  this  method  uses  Dirichlet  conditions  on  the  interior  boundaries,  it 
requires  that  the  subdomains  overlap. 

Van  der  Wijngaart  [7]  revised  the  Schwarz  method  using  asymmetric  in¬ 
terior  boundary  conditions.  If  Neumann  conditions  are  used  for  the  interior 
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boundary  condition  of  Grid  I,  Dirichlet  conditions  are  applied  on  Grid  II.  He 
showed  that  this  treatment  allows  the  requirement  of  subdomain  overlap  to  be 
removed.  More  recently,  Lions  [8]  developed  a  theory  of  the  Schwarz  method 
on  non-overlapping  subdomains. 

A  few  applications  of  this  method  to  fluid  flow  simulation  exist.  Meakin  and 
Street  [9]  simulated  a  3-D  environmental  flow  using  a  composite  grid  method. 
Van  der  Wijngaart  [7]  developed  the  SVVAPR  (Schwarz  Alternating  Procedure 
Revised)  method  which  uses  asymmetric  boundary  conditions  on  the  interior 
boundary. 

In  the  present  work,  we  combine  domain  decomposition  with  the  multigrid 
method.  Henshaw  and  Chesshire  [10]  solved  the  Poisson  equation  using  a  com¬ 
posite  multigrid  technique.  Perng  [11]  and  Perngand  Street  [12]  simulated  flows 
in  complicated  domains  using  a  combination  of  explicit  time  advance  on  individ¬ 
ual  grids  and  a  multigrid  pressure  solver  on  the  complete  composite  grid.  They 
showed  the  effectiveness  of  the  composite  multigrid  method  for  geometrically 
complicated  flows  including  3-D  problems.  However,  they  restricted  the  domain 
decomposition  by  requiring  that  grid  nodes  in  overlapping  domains  belong  to 
both  grids.  This  restriction  reduces  freedom  in  grid  construction  and  is  rmoved 
in  the  present  work.  On  the  interior  boundaries,  Perng  and  Street  adopted  Neu¬ 
mann  boundary  conditions  but  also  used  Dirichlet  conditions. 

Before  describing  the  Navier-Stokes  solver,  we  consider  the  composite  multi¬ 
grid  method  for  a  simple  1-D  model  problem. 

4.1  Multigrid  composite  grid  technique  (1-D  model  problem) 

Consider  the  following  1-D  model  problem. 

y"  =  2,  (0  <  x  <  1),  with  y(0)  =  0,  y'(l)=l  (24) 

The  exact  solution  is 


y  =  x(x  —  1)  (25) 

We  choose  two  overlapping  grids,  one  from  x  =  0  to  x  =  0.6  and  the  other 
covering  x  =  0.4  to  x  =  1;  we  call  the  former,  Grid  I  and  the  latter,  Grid  II. 
Each  grid  is  divided  into  16  equal  elements,  so  the  finest  grid  size  is  0.0375. 
Three  levels  are  used  on  each  grid;  the  finest  has  16  intervals  and  the  coarsest,  4 
intervals.  In  this  1-D  problem,  we  adopt  a  regular  grid,  so  the  coarse  grid  nodes 
coincide  with  the  fine  grid  nodes. 

To  use  Schwarz  iteration,  we  have  to  estimate  the  interior  boundary  value 
by  interpolating  from  the  other  grid.  Two  solution  methods  are  considered.  In 
the  first,  interior  boundary  value  communication  is  limited  to  the  finest  grid 
and  the  coarse  grid  smoothings  are  carried  out  independently.  In  the  second 
method,  data  communication  at  the  interior  boundary  is  allowed  at  every  grid 


level.  Hereafter,  we  call  the  former  the  incomplete  composite  multigrid  (or  sim¬ 
ply  ICMG)  method  and  the  latter  the  complete  composite  multigrid  (or  CCMG) 
method.  Henshaw  et  al.  [10]  used  the  CCMG  method  while  Pemg  and  Street 
[12]  used  the  ICMG  method. 

The  algorithms  for  two  method  are  as  follows. 

a)  ICMG  method 

i)  Iterate  the  equation  on  Grid  I  using  a  local  V-cycle  with  a  guessed  interior 
boundary  value. 

ii)  The  interior  boundary  value  for  the  finest  level  of  Grid  II  is  obtained  by 
interpolation  on  Grid  I. 

iii)  Iterate  the  equation  on  Grid  II  using  a  local  V-cycle. 

iv)  Find  the  boundary  condition  at  the  finest  level  of  Grid  I  by  interpolation  on 
Grid  II. 

v)  Repeat  (i)  -  (iv)  until  the  solution  converges. 

b)  CCMG  method 

i)  Iterate  the  equation  on  the  finest  level  of  Grid  I  using  a  guessed  interior 
boundary  value. 

ii)  The  interior  boundary  condition  for  the  finest  level  of  Grid  II  found  by  in¬ 
terpolation  on  Grid  I. 

iii)  Iterate  the  equation  on  the  finest  level  on  Grid  II. 

iv)  Interpolate  to  find  the  Grid  I  boundary  condition. 

v)  Repeat  (i)-(iv)  (Schwarz  iteration  on  the  finest  grid  level).  There  is  no  need 
to  iterate  to  convergence.  In  the  test  computation,  we  iterated  2  times. 

vi)  The  residual  on  each  composite  grid  is  restricted  to  the  coarser  grid  and  the 
correction  is  iterated  using  procedures  (i)  to  (v). 

vi)  After  solving  at  the  coarsest  grid  level,  work  back  to  the  finer  grids. 

vii)  Iterate  procedures  i)  to  vi)  (V-cycle)  until  the  solution  converges. 

4.2  Interior  boundary  condition  for  multigrid  composite  grid  tech¬ 
nique 

At  the  interior  boundaries,  we  tried  two  types  of  conditions,  Neumann  and 
Dirich  let. 

First,  consider  the  ICMG  method.  In  this  case,  data  communication  is  lim¬ 
ited  to  the  finest  grid  and  the  interior  boundary  condition  is  found  as  follows. 
For  the  Dirichlet  boundary  condition  case, 

Vrifht  =  ITP(yn)  ,  y",  =  ITP(y')  (26) 


8 


where  ITP  is  an  interpolation  operator  and  y1  means  the  solution  on  Grid  I 
and  is  the  right  boundary  value  on  Grid  I. 

We  can  also  introduce  an  overrelaxation  parameter  to  accelerate  convergence 
as  was  done  by  Tang  [13]. 

yii,Kt  =  w/TP(y")  +  (1  -  u,)(t4,M)ol<1  (27) 

In  a  similar  manner,  we  find  yjjj, 

The  case  of  *\  umann  boundary  conditions  is  similar. 

dyli'k  =  ITP(dytl/dx ) 
or 

dyiuht/d*  =  uIT P(dyll/dx)  +  (1  -  w)(dj/' #/,,/dz)o/<f  (23) 

In  the  ICMG  method,  relaxation  on  each  set  of  coarse  grids  is  carried  out 
independently,  thus  the  interior  boundary  condition  for  the  error  is 

flight  =  0  Dirichlet 

dzlight/dx  =  0  Neumann  (29) 

Next,  consider  the  CCMG  method.  The  boundary  condition  on  the  finest 
grid  level  is  found  as  in  the  ICMG  case.  On  the  coarser  grids,  the  treatment  of 
the  boundary  condition  is  somewhat  more  complicated. 

After  iterating  the  correction  equation  on  the  coarse  grid,  the  result  satisfies 
the  following  relation. 

( 1/new  )r.fh(  =  PTP(y^)  Dirichlet 

{dyjttw  /dx)right  =  ITP(dylnltw/dx)  Neumann 

In  other  words,  instead  of  Eq  (29),  the  boundary  condition  on  the  error  on  the 
interior  boundary  should  be 


(ylu  +  *l)n,ht  =  ITP(y’Jd)  +  ITP(iu)  Dirichlet 
(dyiid/dx  +  dil/dx)ritht  =  IT P{dylJd/ dx)  +  ITP(den /dx)  Neumann  (30) 
Again,  we  can  accelerate  the  boundary  conditions: 

il  =  we"  +(1 

where  e*  =  ITP{yIJd)  +  ITP(en)  —  y^j.  The  Neumann  boundary  condition 
can  be  treated  in  a  similar  manner. 

Although  Eq  (30)  is  the  exact  boundary  condition  for  the  error,  we  may 
approximate  it  in  a  way  that  makes  data  communication  at  each  grid  level  in¬ 
dependent,  namely 

elri'ht  =  ITP(in)  or  d'lri,ht/dx  =  ITP(din /dx).  (31) 
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In  the  present  study,  we  tested  both  Eq  (30)  and  Eq  (31). 

4.3  Remarks  on  the  1-D  model  problem 

Next,  we  describe  some  1-D  test  computations.  The  relaxation  method  was 
Gauss-Seidel.  Two  sweeps  were  performed  at  each  grid  level.  On  the  coarsest 
grid,  iteration  was  carried  to  convergence.  For  the  ICMG  method,  one  V-cycle 
was  made  on  each  subgrid  and  the  total  number  of  Schwarz  iterations  was  lim¬ 
ited  to  20.  For  the  CCMG  method,  two  Schwarz  iterations  were  performed  at 
each  grid  level  and  ten  V-cycles  were  carried  out.  Thus,  the  work  is  almost  same 
for  the  two  methods.  B^ause  the  CCMG  method  requires  data  communication 
on  each  grid  level,  its  cost  is  slightly  larger. 

Now  we  present  results  for  the  1-D  test  problem.  The  t^st  computations  are 
(a)  the  ICMG  method  with  Dirichlet  interior  boundary  conditions  (i.b.c.),  (b) 
the  ICMG  method  with  Neumann  i.b.c.,  (c)  the  CCMG  method  with  Dirichlet 
i.b.c.  (Eq  (30)),  (d)  the  CCMG  method  with  Dirichlet  i.b.c.  (Eq  (31)),  (e)  the 
CCMG  method,  with  Neumann  i.b.c.  (Eq  (30))  and  (f)  the  CCMG  method  with 
Neumann  i.b.c.  (Eq  (31)).  We  show  the  normalized  errors  after  20  Schwarz  (1 
V-cycle)  iterations  for  the  ICMG  method  and  after  2  Schwarz  (10  V-cycle)  it¬ 
erations  for  the  CCMG  method  to  judge  the  rapidity  of  convergence.  The  error 
is  defined  by 


tterieal  ye*aet)^ 

(  ’ 

The  errors  are  given  as  functions  of  the  acceleration  parameter  w.  Fig.  3 
shows  the  errors  with  the  Dirichlet  i.b.c.  for  the  ICMG  and  CCMG  methods 
and  Fig.  4  gives  the  errors  with  Neumann  i.b.c.  for  the  same  methods.  Error  I 
is  the  error  on  the  finest  level  of  subgrid  I.  In  both  figures,  results  obtained  with 
both  Eqs  (30)  and  (31)  are  shown  for  the  CCMG  method.  The  ICMG  method 
converges  faster  with  Neumann  i.b.c.  than  with  Dirichlet  i.b.c..  In  both  cases, 
the  convergence  is  accelerated  by  overrelaxation  of  the  interior  boundary  values. 
The  optimal  relaxation  parameters  are  u>  ~  1.6  in  the  Dirichlet  i.b.c.  case  and 
h>  ~  1.3  in  the  Neumann  i.b.c.  case. 

The  ICMG  method  converges  faster  than  the  CCMG  method,  especially 
with  Dirichlet  i.b.cj.  For  both  types  of  i.b.c  j,  Eq  (30)  gives  faster  convergence 
than  does  Eq  (31);  the  difference  in  CPU  time  per  time  step  is  very  slight  so 
the  method  based  on  Eq  (30)  is  faster  overall.  Interior  boundary  value  accel¬ 
eration  is  less  effective  in  the  CCMG  method.  With  Neumann  i.b.c.,  the  best 
convergence  in  the  CCMG  method  is  obtained  with  w  =  1. 

However,  these  trends  change  with  the  external  boundary  conditions  at 
y  =  0  and  1.  For  instance,  wh~n  Dirichlet  boundary  conditions  are  given, 
i.e.,  y(0)  =  y(l)  =  0,  with  Dirichlet  i.b.ca,  the  error  as  a  function  of  u  is  shown 
in  Fig.  5.  The  CCMG  method  converges  much  faster  than  in  the  previous 
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case.  However,  for  optimal  w  (in  this  case  u  =  1.5),  the  ICMG  method  remains 
faster. 

The  convergence  rate  also  depends  on  the  size  of  the  overlap  region.  For 
these  tests,  the  external  boundary  conditions  were  y(0)  =  0,  y'(l)  =  1  and 

Neumann  i.b.c.  were  used.  Teit  computations  were  made  with  20%,  30%, 
40%  and  50%  overlap.  Fig.  6  shows  the  error  for  the  ICMG  method  after  20 
Schwarz  iterations  as  a  function  of  the  overlap.  The  number  of  grid  points  and 
the  number  of  levels  are  the  same  as  before,  16  equal  elements  and  3  multi- 
grid  levels.  The  optimal  acceleration  parameter  is  only  weakly  affected  by  the 
amount  of  overlap,  but  the  convergence  is  much  faster  when  the  overlap  is  large. 

Finany,  we  compared  the  convergence  speed  between  a  single  (non-composite) 
grid  and  a  composite  grid.  We  solved  Eq  (24)  using  32-point  single  grid  and  a 
50%  overlap  composite  grid.  The  finest  level  of  the  composite  grid  has  24  inter¬ 
vals  of  the  same  size  as  the  finest  non-composite  grid.  Three  levels  were  used  in 
both  cases.  The  ICMG  method  was  used  on  the  composite  grid  with  u  =  1.2. 
Fig.  7  shows  the  maximum  residual  on  each  subgrid.  Since  we  must  guess  the 
initial  interior  boundary  conditions  for  the  ICMG  method,  the  residual  at  the 
first  iteration  is  much  bigger  than  in  non-composite  grid  method.  This  requires 
the  ICMG  method  to  take  more  iterations.  However,  the  convergence  rate  is 
almost  same  so  the  cost  increase  is  almost  entirely  due  to  the  overlap. 

4.4  Composite  multigrid  method  for  a  2-D  model  problem 

The  method  can  be  applied  to  2-D  boundary  value  problems.  Since  we  in¬ 
tend  to  apply  the  method  to  the  pressure  equation,  the  Poisson  equation  was 
chosen  as  a  test  problem.  In  the  test  computation,  the  domain  consisted  of  two 
squares  shifted  diagonally  by  40%  of  the  diagonal  length.  The  physical  domain 
and  the  mesh  are  shown  in  Fig.  8. 

Neumann  boundary  conditions  were  applied  at  the  entire  boundary  of  the 
computational  domain.  We  introduced  four  singularities  as  the  forcing;  their 
strengths  add  to  zero  for  consistency  with  the  boundary  condition  d<t>/dn  =  0. 

f  V20dS=  [  } dS  —  [  |^d/  =  0  (33) 

Js  Js  Jl  ° n 

For  the  interior  boundary  condition,  we  can  use  Dirichlet  or  Neumann  con¬ 
ditions.  Because  the  fine  grid  nodes  do  not  coincide  with  nodes  of  the  coarse 
grid,  if  we  use  Dirichlet  boundary  conditions  on  a  staggered  grid,  it  is  difficult 
to  give  boundary  conditions  for  the  error.  Also,  because  Neumann  boundary 
conditions  are  applied  on  the  physical  boundary,  the  programming  is  simpler  if 
Neumann  i.b.c.s  are  used.  Finally,  as  shown  in  the  1-D  test  problem,  Neumann 
i.b.c.s  give  faster  convergence  so  this  is  the  choice  we  shall  make. 

Now  let  us  describe  the  implementation  of  the  interior  Neumann  boundary 
conditions.  As  mentioned  in  Section  2,  we  used  fictitious  points  outside  the 
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computational  domain.  In  Fig.  9,  suppose  points  A  and  B  are  members  of  Grid 
I  contained  in  Grid  II.  In  the  ICMG  method,  values  at  A  and  B  are  obtained 
by  interpolating  Grid  II  data.  The  boundary  condition  for  the  error  requires  its 
gradient  to  be  zero. 

In  the  CCMG  method,  we  need  coarse  grid  boundary  information  as  well. 
To  find  the  boundary  condition  at  point  E  in  Fig.  9,  the  error  at  points  F  and 
G  is  interpolated.  Further,  the  solution  on  the  fine  grid  at  points  A,B,C  and 
D  is  interpolated  to  give  the  derivative  at  points  H  and  I. 

For  interpolation,  we  used  the  Coons  patch  technique  [3]  with  second  order 
accuracy.  Although  analysis  indicates  that  a  third  order  method  should  be  used, 
we  find  that  this  method  is  sufficiently  accurate;  it  does,  however,  lead  to  some 
minor  discrepancies  in  the  results  that  will  be  discussed  below.  Fig.  10  shows 
the  stencil  for  second  order  Coons  patch  interpolation.  P  in  Fig.  10  is  the  point 
at  which  interpolated  data  are  required;  a  and  0  are  the  local  coordinates  of  P. 
The  interpolated  data  can  be  expressed  as 

<j>(P)  =  (1  -  a)D  +  aB  +  (1  -  0)A  +  0C 

-(1  -  a)(l  -  0)4>{a)  -  (1  -  0)a<f>(b)  -  or  0<t>(c)  -  (1  -  a  )04>(d)  (34) 

where  <f>(a )  is  the  value  at  the  corner  a.  The  local  coordinates  (a,  0)  of  the  four 
corners  are  a  =  (0,0),  6  =  (1,0),  c  =  (1,1)  and  d  =  (0,1).  A,B,C  and  D  are 
second  order  polynomials  along  the  edges  ab,  be,  cd  and  da,  respectively  (see  Fig. 
10). 

We  now  show  the  results  of  a  test  computation.  In  Fig.  11,  the  converged 
solution  is  shown.  The  solutions  completely  coincide  in  the  overlap  domain.  In 
Fig.  12,  the  maximum  residual  in  each  subgrid  is  shown  as  a  function  of  itera¬ 
tion  number.  The  fastest  convergence  is  obtained  with  the  ICMG  method  with 
w  —  1.1.  As  in  the  1-D  problem,  the  effect  of  w  on  convergence  is  much  weaker 
for  the  CCMG  method  than  for  the  ICMG  method.  Of  the  interior  boundary 
conditions  applied  to  the  CCMG  method,  Eq.  (30)  yields  the  fastest  conver¬ 
gence.  The  ICMG  method  converges  faster  than  the  CCMG  method.  Thus  we 
used  the  ICMG  method  in  the  Navier-  Stokes  solver. 

To  investigate  the  speed  of  the  composite  grid  method  relative  to  a  single 
grid  method,  we  solved  the  four-singularity  problem  in  i  square  domain.  10 
V-cycle  iterations  were  performed  for  the  single  grid  method  and  10  Schwarz  1 
V-cyde  iterations  were  carried  out  for  the  ICMG  method.  Fig.  13(a)  shows  the 
result  obtained  using  a  single  (33  x  33)  grid  and  Fig.  13(b)  shows  the  result 
using  the  ICMG  method.  The  composite  grids  sue  each  33  x  25  and,  in  the 
overlap  region,  two  subgrids  coincide.  The  two  solutions  agree.  The  single  grid 
method  required  29.1sec  on  a  SUN-3  with  a  floating  point  accelerator  and  the 
ICMG  method  required  61.5sec.  Computations  were  performed  using  double 
precision.  The  maximum  residual  after  10  V-cycle  computations  on  the  single 
grid  is  4.38  x  10-8.  For  the  ICMG  method,  the  maximum  residuals  on  the  two 
grids  after  10  Schwarz  iterations  are  9.43  x  10-6  and  4.20  x  10-6  and  maximum 
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difference  of  the  gradient  on  the  interior  boundary  is  1.17  x  10  8. 

5.  COMPOSITE  MULTIGRID  METHOD  FOR  THE  N-S  EQUA¬ 
TIONS 

The  composite  multigrid  method  is  used  to  solve  the  pressure  equation.  The 
momentum  equations  are  updated  explicitly  on  the  finest  grid.  Because  we  use 
the  QUICK  scheme,  we  need  two  fictitious  points  outside  the  boundaries;  the 
velocity  components  are  obtained  by  interpolation  using  Coons  patch  at  these 
points  (see  Fig.  14).  There  is  a  difference  between  the  stencils  for  pressure  and 
velocity  interpolation;  it  is  shown  in  Fig.  15.  For  pressure  interpolation,  the 
stencil  contains  nine  cells,  while  for  velocity  interpolation,  only  one  cell  used. 
Velocity  data  at  the  cell  corners  are  found  by  extrapolating  the  adjacent  velocity 
data. 

6.  CHECK  OF  NUMERICAL  SCHEME  OF  N-S  SOLVER 


In  order  to  check  the  computer  program,  we  computed  lid  driven  cavity  flow 
at  Re=3200.  First,  we  used  a  single  grid.  The  computations  were  carried  out 
on  33  x  33,  41  x  41,  49  x  49  and  65  x  65  grids  with  A t  =  0.02,0.015,0.015 
and  0.01,  respectively.  The  convergence  criterion  required  the  maximum  in¬ 
crement  of  velocity  in  one  time  step  to  be  less  than  1  x  10-6.  A  CYDRA-5 
mini-supercomputer  was  used;  for  the  33  x  33  case,  9500  time  steps  and  160 
minutes  of  CPU  time  were  required  to  reach  steady  state.  Except  in  the  early 
stages  of  the  calculation,  for  which  the  pressure  solver  was  limited  to  5  V-cycles, 
we  required  that  the  maximum  change  in  the  pressure  between  iterations  be  less 
that  1  x  10-8. 

The  convergence  error  can  be  estimated  as  [14] 


c 


n 


4>n+l  _  0n 
Ai  - 1  • 


(35) 


where  Ai  is  the  largest  eigenvalue  (spectral  radius)  of  the  iteration  matrix.  In 
the  Poisson  solver  with  the  Gauss-Seidel  method,  it  can  be  approximated  by 
Aj  ~  t 2/N7,  N  being  the  number  of  grid  points.  For  the  momentum  equations, 
we  can  not  evaluate  Ai  in  advance  but,  if  the  maximum  eigenvalue  is  real,  we 
can  use  the  estimate  [15] 


Aun+1 

Ai  ~  — - - 

Aun 


(36) 


The  results  show  that  Aj  —  1  ~  O(10-3).  Thus,  we  estimate  the  convergence 
error  as  10“3. 

Since  the  present  scheme  has  second  order  accuracy  in  space,  the  solution 
can  be  expanded  in  the  Taylor  series  about  the  exact  solution. 
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4>  —  4>t*act  +  a(Ax)3  +  6(Ax)3  ^ -  (37) 

On  a  sufficiently  fine  mesh  the  error  should  be  proportional  to  (Ax)3.  Fig.  16 
shows  the  velocity  at  several  points  on  the  vertical  centerline  as  a  function  of  the 
grid  size;  the  error  is  proportional  to  (Aar)3  as  expected.  Results  obtained  by 
Ghia  et  al.  [16]  are  also  shown.  We  estimated  the  exact  velocity  by  Richardson 
extrapolating  the  value  to  Ax  =  0;  the  values  obtained  are  a  little  smaller  in 
absolute  value  than  Ghia  et  al.’s  values.  The  difference  is  largest  near  the  driven 
lid,  but  is  less  than  1.5%.  The  reason  for  these  difference  is  unknown. 

Next  we  show  results  of  the  multigrid  composite  grid  method  for  cavity  flow 
at  Re  ss  1000  on  the  33  x  33  grid.  The  flow  domain  was  decomposed  into  two 
subgrids  of  25  x  33  grid  nodes  each,  the  1CMG  method  was  used  to  solve  the 
pressure  equation  with  the  convergence  criterion  mentioned  above.  In  this  case, 
the  two  grids  coincide  in  the  overlap  region;  the  purpose  of  this  case  is  merely  to 
demonstrate  the  composite  grid  method  and  the  method  is  essentially  the  one 
of  Perng  and  Street  [12].  Fig.  17(a)  shows  the  velocity  profile  at  steady  state 
on  the  single  grid,  while  Fig.  17(b)  shows  the  velocity  distribution  found  by  the 
composite  grid  technique.  In  the  single  grid  case,  2590  time  steps  were  required 
to  reach  steady  state  with  At  =  0.02;  60  time  steps  require  1  minute  of  CPU 
time.  For  the  composite  grid  case,  9  time  steps  require  1  minute  of  CPU  time 
and  2562  iterations  were  needed  to  reach  the  steady  state  with  At  —  0.02.  The 
difference  is  due  to  the  composite  grid  having  1.5  times  as  many  points  as  single 
grid  and,  in  the  single  grid  computations,  only  one  V-cycle  was  allowed  per  one 
time  step,  while  in  the  composite  grid  computations,  we  allowed  five  Schwarz 
iterations  per  time  step.  These  two  profiles  agree  and  there  is  no  discrepancy 
in  the  overlap  region.  Thus  the  accuracy  of  the  scheme  is  confirmed. 

7.  APPLICATION  TO  COMPLEX  FLOW  FIELDS 

Now  we  show  some  results  for  geometrically  complex  flows.  First,  we  solved 
a  lid-driven  two-box  cavity  problem.  The  physical  domain  used  is  the  same  as  in 
the  2-D  two-box  test  computations  (see  Fig.  8).  The  edge  length  and  lid  speed 
are  unity.  The  upper  lid  is  driven  toward  the  right  while  lower  lid  is  driven 
toward  the  left.  The  Reynolds  number  is  1000.  Each  square  had  33  x  33  nodes. 
The  time  step  At  was  0.02,  making  the  Courant  number  less  than  0.64  every¬ 
where.  In  the  computation  of  the  pressure,  the  ICMG  method  with  u>  =  1.1  was 
used  and  five  Schwarz  iterations  with  1  V-cycle  were  allowed  at  each  time  step. 
In  1  minute,  5.6  time  steps  are  taken  on  the  CYDRA-5. 

Fig.  18(a)  shows  the  velocity  profile  at  t  =  6;  there  are  two  symmetric  vor¬ 
tices.  At  t  =  9,  the  two  vortices  have  grown  bigger  and  the  centers  closer,  as 
shown  in  Fig.  18(b).  Eventually  the  two  vortices  merge  into  one  large  vortex. 
The  velocity  profile  at  t  =  12,  shown  in  Fig.  18(c),  contains  one  deformed  large 
vortex.  This  vortex  causes  two  large  recirculating  flow  regions  in  the  corners. 
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Fig.  18(d)  shows  the  steady  state  velocity  profile  and  Fig.  18(e)  shows  the 
corresponding  vorticity  distribution.  Fig.  18(f)  shows  the  pressure  distribution; 
the  pressure  contours  on  the  two  grids  differ  slightly.  The  pressure  is  obtained 
by  solving  Eq  (15);  the  surface  integral  of  the  RHS  of  Eq  (15)  on  the  overlap 
domain  must  be  the  same  on  the  two  subgrids,  i.e., 

I  VUdS=  f  VUdS  (38) 

JSaj 

Since  the  scheme  has  second  order  spatial  accuracy,  Eq  (38)  must  be  obeyed 
to  at  least  second  order  and  the  interpolations  on  the  interior  boundaries  must 
have  accuracy  higher  than  second  oruer.  In  the  present  work,  we  used  the  Coons 
patch  method  with  second  order  peripheral  functions  and  this  is  the  probable 
cause  of  the  small  discrepancies  in  the  pressure  in  the  overlap  region.  The  dis¬ 
crepancies  are  of  the  order  of  the  accuracy  of  the  scheme  and  are  therefore  not 
important. 

Next  we  show  results  for  an  annular-box  combined  cavity  flow.  Fig.  19  shows 
the  flow  domain  and  grid.  The  Reynolds  number  is  1000  based  on  the  length  of 
square  edge.  The  outer  lid  of  the  annular  section  was  driven  counterclockwise 
and  the  left  edge  of  the  square  section  moved  vertically  upward.  Figs.  20(a) 
and  (b)  show  the  velocity  and  vorticity  distributions  at  t  =  49.  Two  large  stable 
vortices  occupy  the  flow  domain.  In  the  overlap  domain,  the  velocities  on  the 
two  grids  agree  very  well.  This  shows  that  the  usefulness  of  the  composite  grid 
strategy  for  geometrically  complicated  flows. 

8.  CONCLUSIONS 

To  simulate  unsteady  flows  in  geometrically  complex  domains,  we  discretized 
the  2-D  unsteady  Navier-Stokes  equations  using  a  staggered  grid  in  curvilinear 
coordinates;  the  accuracy  is  second  order. 

The  effectiveness  of  the  composite  multigrid  approach  for  geometrically  com¬ 
plex  flows  was  demonstrated.  It  may  be  used  both  for  flows  in  which  an  accurate 
time-history  is  required  or  for  the  computation  of  steady  state  flows.  We  investi¬ 
gated  the  convergence  of  the  ICMG  and  CCMG  methods  for  a  1-D  problem  and 
showed  that  acceleration  of  the  interior  boundary  values  is  very  effective.  The 
optimal  value  of  the  acceleration  parameter  depends  on  boundary  conditions. 
In  the  2-D  problem,  the  ICMG  method  again  has  better  convergence  than  the 
CCMG  method;  however,  acceleration  of  the  ICMG  method  is  also  effective. 

Finally,  we  simulated  two  complex  cavity  flows  at  Re  =  1000  to  demonstrate 
the  effectiveness  of  the  method  for  solving  the  Navier-Stokes  equations  in  geo¬ 
metrically  complex  domains. 
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